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Heterogeneity within tumour cell populations is associated with an increase in malignancy and 
appears to play an important role during cancer metastasis. Using in silico experiments, we study 
the interplay between collective behaviours and cell motility heterogeneities in a model system. 
Working with tumour spheroids that contain two non-proliferating cell populations of different 
motile properties, we explore the conditions required for maximal invasion into surrounding tissues. 
We show emerging spatial patterns of cellular organisation and invasion which are consistent with in 
vitro and in vivo observations. This demonstrates that mechanical interactions at the cellular level 
are sufficient to account for many of the observed morphologies of invasion and that heterogeneity 
in cell motility can be more important than average mechanical properties in controlling the fate of 
large cell populations. 

PACS numbers: 87.19.xj, 87.17.Aa, 87.18.Gh. 


Metastasis, the process during which cancer cells mi¬ 
grate away from primary tumour and disseminate in 
other organs, accounts for more than 90% of cancer fatal¬ 
ities PQ. At the cellular level, malignant cell behaviour 
has been associated with an accumulation of gene mu¬ 
tations and a malfunction of numerous regulatory sig¬ 
nalling pathways [2.. Moreover, following earlier pioneer¬ 
ing works [3], recent analyses have established that most 
malignant solid tumours are a mosaic of different cell 
populations with different phenotypic traits [4j. Owing 
to the variability in these factors across different patholo¬ 
gies, patients and even tumours, no common mechanis¬ 
tic pathway leading to the induction and progression of 
metastasis has been yet characterized mm- 
In contrast, recent advances in imaging techniques at the 
tissue scale EM], along with histopathological studies on 
patients’ biopsies [9], draw a surprisingly unified picture 
of cancer invasion in terms of cell migration behaviours. 
In most cancers, despite tremendous variations in gene 
expression and biochemical environment, tumour cells 
are able to invade collectively by maintaining intercel¬ 
lular coordination EHEE], generating hence multicellular 
structures such as cell sheets, strands or clusters remain¬ 
ing cohesive and polarised m- Accumulating evidence 
suggests that the genetic and phenotypic heterogeneity 
observed in these structures plays a key role in the inva¬ 
sion of the tumour cells into surrounding tissues. Collec¬ 
tive invasion of otherwise non-invasive tumour cells can 
be driven by invasion-competent cancer cells [H£E3g or 
stromal fibroblasts PUCE] through a “leader / follower” 
invasion mechanism [9j. These patterns of collective be¬ 
haviour suggest that a system-level analysis is needed to 
account for the generic aspects of tumour cell invasion. 
Cooperative behaviours are observed in a number of bi¬ 
ological contexts such as collective animal behaviours 
na. glassy dynamics of epithelial tissues EMI], wound 
healing [22H24] and early developmental processes ESI 
[2Z]. Moreover, a number of physical factors modulating 


collective behaviours in tissues have also been investi¬ 
gated: geometrical confinement of cells [28H30] . mechani¬ 
cal and topographical properties of the extracellular ma¬ 
trix pJH 133] and forces applied by neighbouring tissues 
[34H37] . Physical models inspired by self-propelled parti¬ 
cles [38, 39] have recapitulated some of these experimen¬ 
tal findings and demonstrated that interactions between 
cell populations and their environments can be under¬ 
stood without resorting to subcellular processes gOl E]. 
But most of these studies have focused so far on well con¬ 
trolled homogeneous cell populations and the role of het¬ 
erogeneity in cell populations remains poorly addressed, 
both experimentally and theoretically. 

In this letter, we use the framework introduced in [42] 
to analyse the effect of cell heterogeneity on tumour in¬ 
vasion. The model involves a Cellular Potts algorithm 
[43] including where appropriate a self-propelled term 
to account for the active motion of the cells. Tumours 
are modelled as 2D packed populations of motile cells of 
spheroid shape surrounded by cohesive tissues made of 
passive cells that mechanically resist invasion. Although 
real tumours are far more complex, this approach pro¬ 
vides a simple dynamic framework to analyse in general 
terms the interplay between cell motility and cell-cell in¬ 
teractions in cell invasion. We will first study the ini¬ 
tiation of invasion in the case of homogeneous tumours 
as a benchmark case. Timescales involved in this study 
are short enough to ignore tumour growth and cell pro¬ 
liferation. We then quantify the effect of population het¬ 
erogeneity by introducing a small proportion of highly 
motile cells into the system. 

In the model, tissues are represented on a 2D lattice 
of 1280 by 1280 pixels, where the value i of each pixel 
codes for the identity of the cell on that location. The 
cell volume Vi(t ) is constrained to a target value Vo [L 2 ] 
by a bulk modulus k [E/L 2 ]. The interfacial cell-cell in¬ 
teractions, such as adhesion or cortical tensions, are all 
accounted for by a single parameter, J [E/L], which sets 
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FIG. 1: Collective invasion for homogeneous tumours, a) 
Snapshots of in silico experiments. Top-left quadrant rep¬ 
resents the initial state of the system (green-blue cells are 
non-motile, yellow-red cells are bulk tumour cells). The three 
other quadrants represent the state of the system for different 
values of fib at t = 4000 MCS. b) T& = f(pb)- Each data point 
is evaluated for a given fib as the mean gradient of A b = f(fib) 
between 0 and 1200 MCS. Inset: Total invasion length of bulk 
tumour cells A b in respect with time, for different values of 
/ ib . Error bars are standard deviations. 


the surface energy of cell membrane. The active motion 
of each cell is driven by a motile force f i (t) = fiiUi acting 
on a fixed substrate, where rii (t) defines the cell polar¬ 
ization axis and fii [E/L\ its active motile force magni¬ 
tude, which depends on the cell type. The motile force 
is used to calculate a migration energy function which 
captures the mechanical work generated by each active 
cell: Wi = —fi'Vi [E\, where Vi represents the position of 
the cell centroid. At any given time, volume constraint, 
cell-cell interactions and active migration energy terms 


can be combined in to an overall energy function, E(t): 

E = £ J (1 - 6 mi(kl) ) + (f (U - Vo) 2 + Wi) 

k,k' i 

(i) 

where i represents the cell index and (k,k f ) represent 
pairs of neighbouring pixels. bi(k)i(k') is 1 when both 
pixels belong to the same cell and 0 otherwise. The sys¬ 
tem dynamics results from the iterative minimisation of 
this energy function through the Metropolis Monte-Carlo 
algorithm [42] . Time is here expressed in Monte Carlo 
steps (MCS), where 1 MCS corresponds to an average of 
one iteration per pixel over the whole lattice. The po¬ 
larization vector of a cell, n^(t), also evolves over time 
and is set along the direction of the mean previous dis¬ 
placements of the cell in the time interval [t — r, t] ; r rep¬ 
resents the time-scale at which cell polarity responds to 
external stimuli and controls the persistence length of the 
cell trajectory in the absence of other cells [44] . To high¬ 
light the role of cell migration, all cell populations have 
the same mechanical, interfacial and dynamical proper¬ 
ties (Vo = 400 pixel 2 , k = 1, J = 5, r m 10 MCS, noise 
level T = 2.5 [E]). They only differ in the value taken by 
their motile forces. Earlier work m has studied in detail 
the role of the different model parameters and identified 
two important behavioural transitions in the populations 
dynamics (c/. the Table 1 of supplementary materials): 

(i) A transition from static epithelial behaviour to col¬ 
lective streaming in homogeneous motile populations, at 
a critical value of the motile force parameter fi c « 0.125. 

(ii) A critical force fi s « 0.225 needed for a single cell to 
migrate through a non motile population. 

We first quantify the short time scale dynamics of inva¬ 
sion in homogeneous tumour spheroids. These structures 
consist of bulk tumour cells of motile force fib and pos¬ 
sess an initial radius Rq of the order of 10 cell diameters 
corresponding to a population of around 300 cells. Inva¬ 
sion is measured as a cumulative length A bi/^b, t) which is 
the sum of the distances to the initial tumour boundary 
for all cancer cells i who left the tumour. This averaged 
quantity encompasses both the spatial extent of invasion 
and the number of invading cells. A& varies almost lin¬ 
early with time for the first thousand MCS as shown in 
the inset of Fig. 0b. The degree of malignancy of the 
system, i.e the capacity of bulk tumours cells to invade 
surrounding tissues, can be characterised as a rate of in¬ 
vasion at short time scales, Tb(fib) = A^/i*,, At)/At. The 
evolution of Tb as a function of the cells’ motile force fib 
is represented on figure FigJTJb, for At = 1200 MCS. 
Homogeneous tumours display behaviours ranging from 
no invasion to weak collective protrusions to widespread 
single cell dispersal (see Fig.l.a). Transitions are con¬ 
trolled by the two critical motile forces fi c and fi s pre¬ 
viously defined (see movies 1 and 2 in supplementary 
materials). These qualitative transitions are reflected in 
the invasion rate as shown in Fig.l.b. The regime where 
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He < Hb < Hs is remarkable as invasion, although slow, 
is observed and takes the form of finger-like collective 
patterns ( cf '. Fig[l]a bottom-right quadrant and movie 
3 in supplementary materials), as described in many in 
vivo studies [7H9]. This suggests that the emergence of 
finger-like protrusions is the result of simple physical in¬ 
teractions in a regime where single cells are not strong 
enough to invade the surrounding tissue, but coordinated 
cell groups can do so by collectively generating larger 
forces. This process occurs without any directional cue or 
specific interaction between cells, other than mechanical 
forces between cells. While this might capture a funda¬ 
mental mechanism involved in the early stages of collec¬ 
tive invasion, the rates are however low compared to what 
is qualitatively expected from the experimental literature 
m and would not be sufficient to explain why these are 
associated with highly malignant situations. The sensi¬ 
tivity and convexity of the invasion rate to the motile 
force of the tumour cells suggests that introducing het¬ 
erogeneities in the motile properties of the cancer cells 
could have a significant effect on invasion rates and pat¬ 
terns. 

We now consider a tumour spheroid surrounded by 
a small number Nf of cells with a larger motile force 
(/i/), inspired by fibroblasts or stronger tumour cells 
sub-populations but representing more generally a cer¬ 
tain level of heterogeneity in the system (cf. Fig. [2ja) . 
We analysed the system’s response for Nf = 6, 12 or 24, 
i.e. about 2%, 4% or 8% of the total motile cell popula¬ 
tion, respectively. The invasive behaviour of bulk tumour 
cells in the presence of these particular cells is quanti¬ 
fied as previously by the invasion length A f(/j,b,Hfi Nf,t) 
and the corresponding invasion rate T f(fib, Hfi Nf) = 
Af(t)/At. To quantify more specifically the effect of 
heterogeneity on the invasion efficiency, we extract the 
percentage increase in the invasion rate per stronger cell: 
Q(/j,b,Hf) = (dN f Ff)/Tb (Fig. [5Jb). As expected, the 
presence of stronger cells brings no significant enhance¬ 
ment to the invasion efficiency for fib > Hs, where the 
bulk of the cancer cells are able to invade on their own 
(see movie 4 in supplementary materials). When both 
bulk and stronger tumour cells generate a motile force 
below the threshold for single cell invasion, there is only 
a marginal effect on the invasion (see movie 5 in sup¬ 
plementary materials). In contrast, when the bulk of 
the tumour is in the regime of collective invasion (/x c < 
Hb < Hs )> the addition of a small number of stronger cells 
(/if > Hs) has a dramatic effect. Each of them increases 
significantly the invasion rate (20% < Q < 45%) for a 
restricted /if value range. The increase in invasion rate 
has an optimum for fib ~ 0.150 and Hf ~ 0-350. To un¬ 
derstand the significance of this optimum, we analysed 
the morphology of the tumour at the onset of invasion. 
Fig. |2ja illustrates the typical spatial arrangements be¬ 
tween stronger and bulk tumour cells in protrusive fin¬ 
gers, in which stronger cells are commonly seen initiat- 
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FIG. 2: a) Snapshots of in silico experiments with 24 stronger 
cells. Top-left quadrant represents the initial state of the sys¬ 
tem (green-blue cells are non-motile, yellow-red cells are bulk 
tumour cells, black cells are stronger tumour cells). The three 
other quadrants represent the states of the system for differ¬ 
ent values of fib and /if at t — 4000 MCS. b) Heat map of the 
normalized invasion rate per stronger cell Q. The maximum 
is observed for /ib — 0.150 and Hf — 0.350. Data set is made 
of 58 points averaged each over 12 different simulations. 


ing and leading fingers (see movie 6 in supplementary 
materials). To quantify these morphological features, a 
graph-based method is used to identify the contact net¬ 
work of cells leaving the tumour. Fingers are defined 
as groups connected to the tumour body and more than 
two cells away from the tumour boundary. The tip is 
defined as the cell of the finger which is topologically the 
further away from the original tumour boundary. The 
length of the finger is the shortest path length between 
its tip and the original tumour boundary, measured in 
number of cells. Fig. [3Ja shows the proportion of fingers 
led by a stronger cell, pooling all data from experiments 
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FIG. 3: Morphology of invading structures in heterogeneous 
tumours, a) Heat map of the proportion of fingers led by 
a stronger cell. Measurements are taken on five time points 
between 3500 and 4500 MCS, and averaged for each on 12 
different seeds. Data set is the same as for Fig. [2]b. b) Im¬ 
pact of tumour heterogeneity on invading structures. (Top) 
Increase in number of fingers per stronger cell with respect to 
fif. (Bottom) Increase in average finger length with respect 
to /if. Circled numbers relate to results in Fig.[2]a. 

with 6 , 12 and 24 stronger cells. In the regime of collec¬ 
tive invasion (fi c < /i 5 < /i s ), up to 70% of the fingers 
can be led by stronger cells, far more than what would 
be expected from their density at the tumour boundary 
(about 25-30%) in average. When fi 5 > fi s , fingers be¬ 
come a transient feature in the invasion and are equally 
likely to be led by any cell type. Consequently, the pro¬ 
portion of fingers led by a stronger cell decreases. Sim¬ 
ilarly, when fib < fi Cl bulk cells cannot migrate on their 
own and any invading structure forming would be pulled 
by stronger cells. The fact that the proportion of fin¬ 
ger led by a stronger cell decreases with fif beyond the 
optimum is however more surprising (see marks 5 and 
6 in Fig. [3ja). Why would an increase in their motility 
reduce their ability to lead fingers? To interpret this, we 


plotted on Fig. 0 b the number of fingers per strong cell 
as a function of /if for fi 5 = 0.150. As /if increases, the 
proportion of strong cells that generate fingers increases 
to a plateau value of approximately one in three. This 
plateau is expected since cell orientations are not biased 
in the model and only a finite fraction of them would 
migrate away from the tumour from the beginning of the 
simulation. The finger initiation step is therefore not 
responsible for the invasiveness decay at large fif. How¬ 
ever, the finger length has a maximum for an interme¬ 
diate value of fif ~ 0.350, consistent with the optimum 
observed in Fig. |2jb. Sustaining the growth of a finger 
is therefore distinct from the ability to initiate it. Both 
properties are however complementary to control tumour 
invasivity; by multiplying the maps of finger number per 
stronger cell and relative increase in finger length ( cf. 
Fig.4 in supplementary materials), we recapitulate the 
profile of Fig. [2]a. What makes fingers shorter for large 
fif is the fact that strong cells tend to escape individu¬ 
ally soon after fingers are initiated, leaving behind short 
fingers that are unable to invade further into the tissue 
(cf. movie 6 in supplementary materials). 

Overall, heterogeneities in the migratory properties of 
cells combined with simple principles of collective mo¬ 
tion appear to be sufficient to reproduce a large number 
of experimentally observed invasion morphologies. The 
model shows that strongly invasive behaviours involve a 
cooperative process between different sub-populations of 
tumour cells. Pattern of cell types along invasive fingers 
emerge without any requirement for specific interactions 
between these cell types. Therefore, understanding the 
role of population heterogeneity goes beyond character¬ 
ising the distribution of their biophysical properties; in 
particular, one cannot assess the invasive and metastatic 
potential of a tumour only based on the proportion and 
strength of its most invasive cells 0. However, the local 
tumour environment and dynamic coupling that exists 
between cells within the population do control morpholo¬ 
gies and malignancy. This paper provides an illustration 
of how diverse the emerging patterns of invasion can be 
with only three cell populations. A general approach 
for the physics of heterogeneities in collective systems is 
still lacking, but there is no doubt that it would provide 
an important contribution to understanding a number 
of biological processes both in animal development and 
clinical research. 
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